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Abstract 

We introduce a monotone (degenerate elliptic) discretization of the Monge-Ampere op¬ 
erator, on domains discretized on cartesian grids. The scheme is consistent provided the 
solution hessian condition number is uniformly bounded. Our approach enjoys the simplic¬ 
ity of the Wide Stencil method |FQ11| . but significantly improves its accuracy using ideas 
from discretizations of optimal transport based on power diagrams |AHA98j . We establish 
the global convergence of a damped Newton solver for the discrete system of equations. 
Numerical experiments, in three dimensions, illustrate the scheme efficiency. 


1 Introduction 


We introduce a discretization of the Monge-Ampere operator, on three dimensional cartesian 
grids, which is simultaneously monotone and consistent. Existing consistent schemes, based 
e.g. on Finite Elements |BN12l INeil2j or Finite Differences |LR05] . are not monotone, and 
thus require the PDE solution to be sufficiently smooth, and the numerical solver to be well 
initialized. Existing monotone schemes, based on Wide Stencil discretizations (FOlll i()be06|, 
suffer from a consistency error depending on the discretization stencil angular resolution. Filtered 
schemes |FQ13j combine a monotone and a consistent scheme, and attempt to cumulate their 
robustness and accuracy; improving either of the constituting schemes will benefit to the filtered 
combination. A monotone and consistent scheme is introduced in (BCM14) . but it is limited to 
two dimensions. Geometric approaches |OP89| are discussed in the third paragraph. 

The proposed numerical scheme belongs to the Wide-Stencil category |Qbe06j . in the sense 
that we actually define a family of schemes parameterized by a user chosen stencil. Larger stencils 
provide consistency for strongly anisotropic problems (i.e. for which the solution hessian is almost 
degenerate), but at the cost of an increased computation time. The choice of stencil is left to 
the user; let us mention that in the special case of |BCM14| an automatic (solution adaptive, 
local, anisotropic, and parameter free) stencil construction could be designed. Our numerical 
experiments show that small stencils, of radius \/3 or vA see the table page 13 already yield 
convincing results. The scheme is dimension independent, but we emphasize its application to 
three dimensional domains, which is tractable and tested. 

Our approach is also inspired by |OP89| and the discretizations of optimal transport |AH A981 
IMerlllFLevl4j based on global geometric structures, up to two differences. The first modification, 


a symmetrization see Remark 1.9 is required to operate our method with the Dirichlet boundary 
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conditions of the standard Monge-Ampere problem ([!]), instead of the second boundary condi¬ 
tions implicit in optimal transportation. The second modification localizes these methods by 


limiting interactions to close discretization points, see Remark 1.10, which considerably simpli¬ 
fies their numerical implementation. The methods |OP89l IAHA981 IMerlll ILevl4| indeed rely on 
global geometric structures named power diagrams, which generalize Voronoi diagrams. Their 
construction requires state of the art methods of discrete geometry, which especially in 3D are 
still an active subiect of research. For instance [Lev m mentions arbitrary precision arithmetic, 
arithmetic filtering, expansion arithmetics and symbolic perturbation, merely for the consistent 
evaluation of geometric predicates. Our approach is in contrast local and requires none of these 
subtleties. We show that this simplification preserves consistency in the setting of viscosity 
solutions, see §2.3| but at the following price: the solution hessian condition number must be 
uniformly bounded, and the weak Alexandroff solutions cannot be recovered. 

We fix throughout this paper an open, convex and bounded domain 12 C M rf , d > 2 (d = 3 
in the numerical section ^3|. Given a positive density p G C 0 (12,Mlj_), and some Dirichlet data 
a g c°(an,R), we set the goal of approximating numerically the unique viscosity solution 
|GTL92I iGMOTj of 

det(V 2 «) = p on 12, 

u = a on <912, (1) 

u convex, 

where V 2 u denotes the hessian matrix of u. The PDE domain 12 is discretized on a cartesian grid 
X. Up to a linear change of coordinates, encoding scaling, rotation and offset, we may assume 
that 

icfinz' 1 . 


The discussion of the discretization of the boundary <912 is postponed to Remark 2.4 

Definition 1 . 1 . We denote by HJ the collection of maps u : X U <912 —> M. 

Definition 1 . 2 . A discrete operator D associates to each u G U a discrete map Du : X —> M. 

Given some discretization D of the Monge-Ampere operator, the counterpart of ([I]) takes the 
form: find u G U such that 

{Vu = p on X, 

|'fi = a on <912. 

The constraint “u convex” is not spelled explicitly in ([2]), contrary to ([!]), but some discrete 
counterpart of it often follows from the identity Du = p |BCM14| . Before discretizing the 
Monge-Ampere operator, we need to introduce a more basic operator A e , e G Z d , aimed at 
approximating the second order difference (e, (V 2 rt(x))e). If s G T, x + e G X and x — e G X, 
then we set classically 

A e u{x) := u(x + e) — 2 u(x) + u(x — e). (3) 

In general, one may not have x + e € X, for instance if x + e lies outside 12. Hence we introduce 


h e x := min{/t > 0; x + he G X U <912}. 


(4) 


Let h + := h x and h := h x e . Only one linear combination of u(x), u(x + h + e) and u(x — he) 
is consistent with (e, (V 2 u(x))e), namely 


A e u(x) := 


h+ + h~ 


u(x + h + e) — u(x) u(x — h e) — u(x) 
h+ + h 1 


(5) 
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We next illustrate the notions of consistency and monotonicity (degenerate ellipticity) using 
two now classical discretizations of the Monge-Ampere operator. The Finite Differences (FD) 
[ LR~05j approximation of the Monge-Ampere operator is defined as the determinant of a naive but 
consistent approximation of the hessian matrix of u by finite differences: denoting by (ej)i<j<d 
the canonical basis of 


V FD u(x) := det(<%)i<jj< 3 , with d tJ 


A ei u(x) if i = j, 

i (A ei+ej u(x) - A ei _ ej u(x)) if i / j. 


( 6 ) 


The Wide Stencil (WS) fFOTT] approximation of the Monge-Ampere operator is defined as 
follows: denoting by B C (fL d ) d a hnite collection of d-plets of pairwise orthogonal vectors 


TXS s u(x) := min TT 

B&B 11 


max{A e u(x), 0} 


ess 


(7) 


Definition 1.3. An operator T> is said Degenerate Elliptic of second order (DE2), with stencil V, 
iff for all x G X, Du{x) is a non-decreasing function of the second order differences (A e u(x)) eg y. 

By construction, scheme WS is DE2, with stencil {e G Z d ; 3B G B, e G B}. Scheme FD in 
contrast is not DE2. Degenerate Ellipticity comes with strong guarantees including a maximum 
principle for the discrete solutions of ([2]), and guaranteed convergence of Euler iterative solvers 
for this discrete system |Qbe06| . It allows to recover non-smooth viscosity solutions of ([!]), see 
, and plays a crucial role in the proof §2.2| of the global convergence of a damped Newton 
solver for (j2|. 

Let Sd denote the set of symmetric dx d matrices, and let S'j" C Sd be the subset of positive 
definite matrices. In order to analyze the consistency of these schemes, we introduce for each 
M G S'j" the norm || ■ \\m and the map um G U defined by 


£ 2.3 


|*||M ■= \J(x,Mx), 


um{x) 



( 8 ) 


Definition 1.4. The consistency set of an operator T> is the collection of all M G such that 
T>um = det(M), identically on X. 

The consistency set of scheme FD is the whole ; in fact, the identity Dum(x) = det(M) also 
holds for non-positive matrices, although they are irrelevant for our application. The consistency 
set of scheme WS is in contrast of empty interior ; precisely, it consists of matrices M G S'j - for 
which some B G B is an eigenbasis [Qbe06| . Stated otherwise, the consistency of scheme WS is 
an asymptotic property, obtained by growing the stencil size to infinity. 

A combination Vu(x) = w(u,x)Vg S u(x) + (1 — w[u,x))V FY> u(x) of schemes FD and WS is 
considered in |FQ13j . The weight w(u, x) G [0,1] depends on the local behavior of u close to x, as 
well as on the discretization scale and the angular resolution of the stencil B. Strictly speaking, 
the weighted scheme loses both the degenerate ellipticity of Dff s , and the consistency of D FD . 
Nevertheless |FQ13| establishes convergence in the setting of viscosity solutions, and illustrates 
numerically (in two space dimensions) that the weight w(u,x ) favors the consistent scheme V FD 
except close to the most singular features of the solution. Second order accuracy may in principle 
be achieved even on a degenerate equation, contrary to the method proposed below. The filtered 
scheme construction proposed in |FQ13| is quite flexible, and could be improved by replacing 
the monotone scheme WS with the more accurate, and still monotone, scheme proposed in this 
paper. 
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We propose an alternative solution to the apparent conflict between between degenerate 
ellipticity and consistency. The coordinates of a vector e = (ai, • • • ,a d ) G are said co-prime 
iff gcd(ai, ■■■ ,a d ) = 1. 

Definition 1.5. We limit our attentions to stencils F C Z d \ {0} which are finite, symmetric 
w.r.t the origin, span M. d , and which elements have co-prime coordinates. The proposed operator 
is 

Vyu(x) := Lebjg G M® Ve G V, 2(g,e) < A e u(x)}, (9) 

where Leb denotes d-dimensional Lebesgue measure. 


The Degenerate Ellipticity of T>y is clear: if any second order difference A e u(x) increases, 
then the convex polytope appearing in Q increases for inclusion, hence also in volume. Refining 
slightly this argument we obtain the derivative of T>yu(x ) with respect to A e u{x): namely the 
(d — l)-dimensional measure of the facets of (J9| defined by the equality constraint 2\(g,e)\ = 
Ae'u(.x), divided by ||e|| (as a result, Vyu(x ) is continuously differentiable in u, in contrast with 
0 )- Computing polytopes defined by linear inequalities like ([9| is, by convex duality, equivalent 
to computing the convex envelope of a set of points |PS12| . Numerous computer libraries are 
available for that purpose, such as TetGen®or CGAL®(the author made his own routine which 
takes advantage of the symmetry of the polytope ([9])). 

The expression ([9j seems more complex to evaluate than (|6]) or ([7]). Yet in our numerical 
experiments ^3j the cost of evaluating the operator T>y was dominated by the cost of solving linear 
systems (when solving (J2]) with a damped Newton solver), for which the MUMPS®library was 
used. The applicability of Newton’s method to the Monge-Ampere problem is also investigated 
in [LR.05]. for the continuous problem, and in [FOllJ and [Nei m for some discretizations. Its 
super-linear rate of convergence, in the neighborhood of the problem solution, makes it appealing 
for applications. In contrast, gradient descent |AHA98| or first order Euler |Qbe06] schemes 
converge slower but benefit from global convergence guarantees with monotone discretizations 
(iterates converge to the unique solution of the discrete problem, independently of initialization). 


For the discretization of interest, we prove i2.2 that a damped Newton method benefits from 
both a super-exponential rate of convergence close to the discrete problem solution, and a global 
convergence guarantee. 

The following results, established §2.1| show that the consistency set of the proposed scheme 
is of non-empty interior, in contrast with scheme WS. By Corollary |1.8[ a finite stencil is sufficient 
to guarantee consistency for all matrices M G Sj - with condition number below a given bound. 


Definition 1.6. For each matrix M G Si~, we introduce the Voronoi cell and facet 


Vor(M) := {g G M® Ve G Z d , || 5 || M < h ~ e\\ M }, 
Vor(M; e) := {g G Vor(M); \\g\\ M = ||9 - e\\ M }- 


An M-Voronoi vector is an element e G Z rf \ {0} such that Vor (M; e) 0; it is said strict iff the 

facet Vor(M;e) is (d — 1)-dimensional. 


Proposition 1.7 (Consistency). A matrix M is in the consistency set ofDy iff V contains all 
strict M- Voronoi vectors. 

Corollary 1.8 (Finite stencils are enough). Let k > 1 and let V collect all elements e G Tfi 
with norm ||e|| < ny/d and co-prime coordinates. Then the consistency set of Vy contains all 
M G such that k > y / ||M||||M _1 ||. 
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Figure 1: Voronoi cell Vor(M), Voronoi vectors e, and ellipse {x £ M d ; (x,Mx) < 1}. Left: 2D, 
Right: 3D. Note that e traverses the facet Vor(M;e) at the point e/2. 


Outline: The results on consistency are established §2.3| and the global convergence of the 

damped Newton solver in (2.2 We also study in §2.3| convergence as the grid scale tends to zero, 
in the setting of viscosity solutions. Numerical experiments ([3] illustrate the method’s efficiency. 

Remark 1.9 (Symmetrization). A non-symmetrical variant of our Monge-Ampere operator dis¬ 
cretization (|9|) can be defined as follows: for x £ X such that x + e £ X for all e £V 


'Dyu(x) := Leb {<7 £ M d ; Ve £ V, (g , e) < u(x + e) — ii(:r)}. 


( 10 ) 


In the case of a quadratic function um, M £ Sf, the polytope appearing in (10) is merely a 


translation of ®> by the vector Mx, so that VyiiM = T>yUM ■ dn the case of a general u, the 
polytope (10) can be non-symmetric, in contrast with ©• This asymmetry can actually help 
extract weak Alexandroff solutions of ([!]) with non-symmetric subgradient sets, such as the two 
Diracs test case in wm . 

The drawback of the expression © is that it is only correctly defined sufficiently far from 
the boundary dLl. There is no simple way, to the knowledge of the author, to incorporate in (10) 
the values of u on <9f l, and remain consistent with the Monge-Ampere operator det(V 2 u) in the 


strong sense of Definition l.f. Hence our choice of defined in terms of the symmetric second 
order differences A e u(x), e £ V. 

Remark 1.10 (Localization). Consider the point dependent, largest possible stencil V(x) = {e £ 
Z d ; x + e £ X}. Then Dy^u(x), see ( |10[ ) ; is the measure of the power cell associated to x in the 
power diagram based discretization ]OP8fA[AHA98[\Merll[ Levif l of optimal transport. Since the 


stencils ( V(x)) X £x ar e typically huge, it is not reasonable to evaluate V'y^u(x) independently 
for each x. Instead a global power diagram needs to be constructed, using complex geometric 
procedures which validity requires extremely careful evaluations of geometric predicates \Lev W- 


2 Proofs of the main results 


We establish in ^2.1 the numerical scheme consistency, show in (2.2 the global convergence of 


a damped Newton solver for the discrete system, and prove in (2.3 a convergence result in the 
setting of viscosity solutions as the discretization grid scale tends to zero. 


2.1 Consistency 

Our consistency analysis relies on elementary properties of Voronoi cells of lattices, see Definition 


1.6 


and [CS92 J for more on this topic. Our first step is to bound their diameter. Let n(M) 2 
denote the condition number of a matrix M £ Sjf , i.e. k(M) := v/pff jM^l 
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Lemma 2.1. Let M E S£. Any point g E Vor(M) satisfies ||g|| < \n(M)\fd. Any M-Voronoi 
vector e satisfies ||e|| < and has co-prime coordinates. 

Proof. Bound on the norm. Let g E Vor(M), and let e g be obtained by rounding the coordinates 
of g to the nearest integer, so that || g — e g \\ < \\[d. One has A~||g|| < \\g\\M < 115 — e g \\M < 
X + \\g — e g \\, where A _ ,A + are respectively the smallest and largest eigenvalues of Ms. Thus 
IMI < {\ + /\~)\\fd as announced. If e is an M-Voronoi vector, then there exists g E Vor(M) 
such that ||p|| = \\e — g ||, hence ||e|| < 2||g|| < n(M)Vd. 

Coordinates are co-prime. Consider a vector which coordinates are not co-prime: ke, with 
k > 2 and e E Z d \ {0}. Then for any g E one has || ke — g\\ 2 M + (k — 1)H^H^f = ^’ll e — 
g \\ 2 m + (k 2 — fc)||e|||^. Hence ||e — g\\M < max{||fce — g\\M, HffHiVf} and therefore ke cannot be an 
M-Voronoi vector. □ 


Corollary 2.2. Let M E S'V, and let E be the collection of strict M-Voronoi vectors. For any 
set V C 7L d one has Vor(M) C {g E M d ; Ve E V, 2 (g,Me) < ||e|||^} ; with equality iff E Cf. 


2 

Mi 


Proof. For any g, e E M , we have the equivalence \\g\\M < ||5 — e|| m 2 (g,Me) < ||e 

obtained by squaring both sides of the first inequality. Hence Vor(M) is a convex polytope, 
defined by a family of linear inequalities indexed by e E By Lemma 


2.1 


one can eliminate 

all inequalities but a finite number. Among these inequalities, only those corresponding to strict 
M-Voronoi vectors are active, in the sense that they define a facet of Vor(M). The result 
follows. □ 


We next identify the volume of a Voronoi cell, and use it to establish our scheme consistency. 


Lemma 2.3. For any M E one has Leb(Vor(M)) = 1. 

Proof. The set Vor(M) collects elements g E which are closer to the origin than to any other 
point e E Z d . As a result the translates of Vor(M), by all offsets e E Z d , tile the space W l up to a 
negligible set. Thus Leb(Vor(M)) is the co-volume of the lattice Z d , namely 1 as announced. □ 


Proof of Proposition 1.1 (Consistency). Let M E S([ , let x E X, and let um be the quadratic 
map d 8 l). By construction A e UM(x) = ||e|||^, for any e E Z d . Hence Vyum{x) is the volume of: 


{g E M rf ; Ve E V, 2(g,e) < ||e||^} — M{g E 


*; Ve E V, 2( 5 , Me) < ||e||^} D MVor(M). 


where we used Corollary 2.2 and abusively denoted by ME := {Me; e E E} the image of a set 


E by the linear map M. Therefore Vvum(x) > det(M) Leb(Vor(M)) = det(M) by Lemma 2.3 


Equality holds iff the above inclusion of polytopes is an equality, equivalently iff V contains all 
strict M-Voronoi vectors by Corollary |2.2| □ 


Proof of Corollary 1.8 (Finite stencils are enough). Combine Proposition 1.7 and Lemma 2.1 


n 


2.2 Global convergence of a damped Newton solver 

We establish the convergence of a damped Newton solver for the discrete system (J2|. Before 
doing so, we need to clarify the implementation of boundary conditions. The stencil V is fixed 


in the following, and obeys the properties of Definition 1.5 
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Remark 2.4 (Boundary discretization). Maps m£D are in principle, see Definition ]! . 1\ defined 
both on the discrete sampling X of Cl, and on the whole uncountable boundary dCl. Fortunately, 
only finitely many values of u on dCl play an active role in the system of equations Q: those on 

dX := {x + h%e\ x G X, e G V} n dCl. 


This discretization, with ff(X U dX) unkowns, is used in our experiments, and in the proof of 
the damped Newton solver global convergence below. 

Alternatively one may (i) not introduce any unknown for the boundary values, but use the 
given Dirichlet data a as in IBCMlfitf , or (ii) introduce some unknowns on an arbitrary boundary 
sampling dX', extended to dCl\dX' by some interpolation procedure. Unfortunately (i) leads to 
initialization difficulties for the iterative solver, since one needs to find a strictly convex seed uq 
with prescribed boundary values on dCl, and (ii) may limit the accuracy of the scheme if first order 
interpolation is used, or violate its degenerate ellipticity in the case of high order interpolation. 

Our first step is to establish, in Corollary |2.6[ that the Jacobian matrix associated to the 
discrete Monge-Ampere system ([2j is invertible. This follows from the invertibihty of diagonally 
dominant matrices, recalled in the next lemma, and from the degenerate ellipticity of the proposed 
operator Vy. 

Lemma 2.5. Let A be an n x n matrix such that for each 1 < i < n 

1^*1 >£ I Ail- (11) 

j A* 

Assume that for each index 1 < io < n there exists a chain i$, ■ ■ ■ , if, such that i^ satisfies strictly 
inequality and Ai rt i r+1 0 for all 0 < r < k. Then A is invertible. 

Proof. For each index 1 < i < n, let k(i) denote the length of the smallest chain of indices as 
above. Let x G be such that Ax = 0. Among the indices 1 < i < n such that the vector 
component |xj| is maximal, choose one which minimizes k(i). From (Ax)i = 0 we obtain 


Au\\xi\ < ^2 \ A ij\\ x j\> tllus (\ A ii\ - ^2 1 AiDI^I + l^-KN - Nil) = U + V < 0. 

3 A* 3& 3Ai 


Both terms U and V are non-negative by construction, hence U = V = 0. If k(i) = 0, then from 
U = 0 we obtain \xf\ = 0, hence x = 0. If k{i) 0, then from V = 0 we obtain \xj\ = |xj| for 
all indices 1 < j < n such that A,j 0. One of these indices satisfies k(j) = k(i) — 1, which 
contradicts our choice of i. This concludes the proof. □ 

Let Uq be the collection of maps u : X U dX —> M such that: 


Vx G X, Vyu{x) > 0. (Equivalently: Vx el, Ve G V, A e u(x) > 0.) 
Let / : Uq —> M' Yuax be dehned by 


f(u)(x) 


In Vyu{x) if x G X , 
u(x) if x G dX. 


( 12 ) 


( 13 ) 


Corollary 2.6. For each u G Uq, the Jacobian matrix df(u) is invertible. 
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Proof. Let m £ Do and let A := df(u), which formally is a matrix associating a coefficient to 
each index pair (x, y), x, y E X U dX. The line of A corresponding to each x E dX has a single 
non-zero coefficient, namely 1 at index (x, x ), hence this line satisfies strictly. By degenerate 
ellipticity of XV, the line corresponding to any x E X satisfies ©• We prove in the following 
the chain property of Lemma 2.5 which implies the announced invertibility. 

Let x E X, and let K x be the polytope appearing in ([ 9 ]), which by construction is convex, 
compact and symmetric w.r.t the origin. For each e E V let F e denote the facet of K x defined 
by the equality constraint 2(< 7 , e) = A e u(x). Let V x : = {e E V; A e > 0} where A e denotes the 
(d — l)-dimensional measure of F e . Since the polytope is symmetric, one has A e = A_ e for all 
e £ V, hence V x is symmetric. Since the polytope is compact, and since the exterior normal to 
F e is e/||e||, the set V x spans M d . The coefficient of the Jacobian matrix A at index (x,x + h%e) is 
4Ae/mK + h-*)\\e\\ Leb(K x )), using ([ 5 ]) and the geometric argument that the polytope volume 
variation is at hrst order given by the facets areas A e + A_ e times their normal displacement. 
This coefficient is hence positive if e E V x . 

Among the points x + h x e, e E V x one at least is strictly closer to dLl than x. By induction 
we can thus build a finite chain x = xo, ■ ■ • , Xk-i E X such that Xk E dX and the coefficient of 
A at index (x r ,x r +i) is non-zero for each 0 < r < k. The announced result then follows from 
Lemma 12.51 □ 


Our second step is to show that (13) is a proper map: the preimage of any compact set is a 


compact set, a property which is tightly linked with the well-posedness of the PDE ([!]). Let 

8 X U := {g E M d ; Vy E U, (g, y - x) < U{y) - U(x)} (14) 

denote the subgradient of a convex map U E at a point x E fi. We connect in Lemma 


2.7 the proposed operator T>yu(x) with the Lebesgue measure of d x U, where U is the lower 


convex envelope of u. The properness of (13) then follows in Corollary 2.9 from the maximum 
principle of Alexandroff-Bakelman-Pucci. 

Lemma 2.7. Let u E Uo, and let U : Hull(X U dX) —> M be the maximal convex map bounded 
above by u. Then Leb(<9 x ?7) < h^Vyu(x) for all x E X, with h* := max{/&®; x E X, e E V}. 

Proof. By construction U < u, and for any x E X such that U(x) < u(x) one has Leb (d x U) = 0, 
so that the announced inequality holds. Assume that U(x) = u(x), and let g + ,g~ E d x U. Let 
e E V and let h + := h%, h~ := h~ e as in ([ 5 ]) . One has ( g + ,h + e) < u(x + h + e) — u(x), and 
(g ~, —h~e) < u(x — h~e) — u(x) by (14), thus 


(g + -g-,e) < 


u(x + h + e) — u(x) u{x — he) — u(x) h + + h 


h+ 


+ 


h~ 


A e u(x) < h*A e u(x) 


Therefore \{d x U — d x U) C h*K x , where the left hand side is a Minkowski sum of sets, and K x 
is the polytope appearing in ([9|. The announced estimate then follows from this inclusion and 
Brunn-Minkowski’s inequality. Q’ 

Theorem 2.8 (Alexandroff-Bakelman-Pucci’s maximum principle [Gut01] ). Let u E and 

let U be the maximal convex map bounded above by u. Then 


min u — rrhn u < (Leb(G) /bJd) d diam(fl), 

dn n 


with G := d x U, 


(15) 


u(x)=U ( x ) 


where diarn(Q) := max{||x — y||; x, y E fl}, and cod is the volume of the d-dimensional unit ball. 
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Corollary 2.9. The map f : Uq —> M Au9A " is proper. 


Proof. Since / is continuous, it suffices to bound by above and below the values of an arbitrary 
m 6 Do in terms of those of f(u). We have the obvious estimate f(u)(x) = u{x ) for values at 
x G dX. 

Lower bound on X. Let Q' be the interior of Hull(X U dX), and let u' G C°(n ,M) be 
defined by u'(z) = min^gxuax u(x) + A||x — z\\, where A > 0. Let U G C°(n , M) be the 
maximal convex function bounded above by v!. If the constant A is sufficiently large, then U 
also is the maximal convex function bounded above by u. In addition for all s 6 fl', one has 
U(x) = u'{x) x G X U dX. We obtain applying Theorem 2.8 the desired lower bound on 
min{u(x); x G X}, since 


min u' 
Q’ 


min u, mini/ = min u, Leb(G) = Leb(d x U) < hi Vyuix) 
xudx ’an' ox v ^ v x > - * v ' 

x&X x£X 


hi Y e f(u)(x) - 

xEX 


Upper bound on X. Let x G X U dX be such that u(x) is maximal, and let us assume for 
contradiction that x ^ dX. Since A e u(x) > 0, for any e G V, one has either u(x + h + e) > u(x), 
or u(x — h~e) > u(x), with the notations of (|5]). This contradicts the maximality of u(x), and 
concludes the proof. □ 


Numerous variants of the Newton method exist |EW94| : an elementary one is presented be¬ 
low for completeness, which guarantees global convergence for the system (|2]) of interest. This 
damped Newton algorithm is an iterative equation solver, which recursion rule © involves an 
adaptively chosen “damping” parameter 5. In practice 6 is typically small in the first iterations, 
and equal to 1 in the last ones, which coincide with the classical Newton method and enjoy 
quadratic convergence. Convergence is guaranteed for maps which, like ([5]), are shown to be 
proper and at each point a local diffeomorphism. Note that these assumptions, plus the con¬ 
nectedness of the source domain and the simple connectedness of the target domain (here both 
satisfied), imply by Hadamard-Levy’s theorem that / is a global diffeomorphism. 

Proposition 2.10 (Global convergence of the damped Newton algorithm). Let N > 0, let 
Uo C M. n be an open set, and let f G C' 1 (Uo, M. N ). Assume that f is proper and that the Jacobian 
matrix F(x) := df(x) is invertible for each x G Uo- Consider y G M. N , and for each x G Uo, 
6 G [0,1], the Newton update 

Newton(x, S) := x + 5F(x)^ 1 (y — f(x)). (16) 

Let xo G Uo, and for each n > 0 let x n +i := Newton(x ra , S n ), where 5 n = 2~ kn and k n is the 
smallest non-negative integer such that \\y — f{x n+ i)\\ < (1 — S n /2)\\y — f(x n )\\ . Then f (x n ) y 
as n —> oo. 


Proof. Introduce the set K := {x G Uo; || y — /(x)|| < \\y — /(xo)||}, which is compact since / is 
proper. For any x G Uo one has the Taylor expansion: for small 5 

/(Newton(x, 5)) = f(x) + 5F(x)F(x)~ 1 (y - f(x)) + o(5||F(x) _1 (?/ - f(x)) ||) 

= V ~ (1 - 8)(y - f(x)) + o(5\\y - f(x )||). (17) 

Hence ||y — /(Newton(x, 5))|| = (1 — 6 + o(S))\\y — /(x)|| using © , which is smaller than 
(1 — 6/2)\\y — /(x)|| for an open range of 6 G]0, A(x)[. By compactness, this property holds for 
all iGif with an uniform open range ]0, A[. 

Thus So is well defined, bounded below by A/2, and by construction \\y — /(xi)|| < (1 — 
A/2)||y — f{x o)|| so that x\ G K. The result follows from an immediate induction argument. □ 
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2.3 Convergence in the setting of viscosity solutions 

In this section, we let the discretization grid scale tend to zero, and study the convergence of the 
minimizers to the discrete problems (|2]). For all integers n > 1, let X n := flfln -1 Z d . Let HJ n be 

the collection of (semi-)discrete maps u : X n U dll — > M. The stencil V is fixed and obeys the 

assumptions of Definition |1.5| For any u E HJ n , x E X n , and M E Sd let 

V n (u) := Leb{<? E M d ; Ve E V, 2(g,e) < A e/n u(x)}, 

D(M) := Leb {g E M d ; Ve E V, 2 (g, e) < (e, Me)}. (18) 

We denoted by A e / n u(x) := n 2 (u(x + e/n) — 2u(x) + u(x — e/n)) the standard approximation 
of (e, ( V 2 u(x))e ) on the grid X n (This expression assumes that x ± e/n E X n , and should 
be modified as in ([5]) otherwise). Given a density p E C'°(fl,M^_), and some Dirichlet data 
a E C°(dfi,R), we study the problems 

u : —> M, 

D(X 2 u) = p on D, (19) 

u = 1 7 on dQ. 


u E 

< V n (u ) = p on fl, < 

u = a on dVt. 


2.15 


By 12.2, the discretized problem (19 left) has a unique solution u n E U n . We show in Corollary 

Importantly, if X 2 u oc exists at each 
is also the unique solution to the 


that (19 right) also admits a unique solution u c 


x 


E and belongs to the consistency set of T>y. then u a 
Monge-Ampere equation (JTj). The study of (19 right) relies on the concept of viscosity solutions 
|CTL92j . 


Definition 2.11. A sub-solution (resp. super-solution) of \19\ right) is an upper-semi-continuous 
(resp. lower-semi-continuous) map u : D —> M such that (I) u < a (resp. u > a) on dQ and (II) 
for any iESI and any ip E C 2 (fl) such that u — p has a local maximum (resp. minimum) at x, 
one has D(\7 2 ip(x)) > p(x) (resp. < p(x)). 

A solution to (19 right) is a map u : D —> M which is both a super-solution and a sub-solution. 


Replacing “local maximum” with “strict global maximum” (resp. minimum) in Definition 
2.11| y ields an equivalent notion of sub- and super-solutions [ 011192 ]. Super- and sub-solutions 
of (|19| right) obey a comparison principle, proved Proposition 2.13 Given symmetric matrices 
M, M' E Sd, we write M A M' iff the difference M' — M is positive semi-definite. 

Lemma 2.12. For any M, M' E Sd, one has: M A M' => D(M) < D(M'). More quantitatively, 
if M, H E Sd are such that D(M) > 0 and D(H ) > 0 then 

D(M + H)-d > D{M)-d +D(H)l. 


Proof. Let K(M ) C 

M A M' => Ve E V, ||e||^ A 


be the polytope appearing in (18), so that D(M) = Leb K(M). Then 

D(M) < D(M'). 


| p ||2 


K(M) C K(M') 


Second point: since (e, ( M+H)e) = (e, Me) + (e, He), the set K(M+H ) contains the Minkowski 
sum K(M) + K(H). The announced result follows from Brunn-Minkowski’s inequality. □ 


Proposition 2.13. If u* 


is a sub-solution of (19, right), and u* 


a super-solution, then u * < u*. 
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Proof. The result does not immediately follow from f CIL92j because the operator —D{S7 2 u) is 
only degenerate elliptic. The following operator is in contrast strictly elliptic when e > 0: 


F e (u) : = —D(V 2 u)<i + eu. 


( 20 ) 


Let M := sup^ u* — u*, which is finite by upper-semi-continuity. For contradiction, we assume 
M > 0. Let also r := max{||x||;x G $2} and 


u £ (x) := u*(x) + ^(||.t|| 2 - r 2 ). 


One has u £ < u* < u* on dui, and by Lemma 2.12 


F e (u*) - F £ (u £ ) > eM - e(u* - u*) > 0. 


( 21 ) 


( 22 ) 


Applying the standard comparison principle [CIL92] to the strictly elliptic F e we obtain u* > u £ , 
hence u* — u* > —eMr 2 / 2. Letting e — y 0 yields u* > u* as announced. □ 

Weal solutions can be extracted from sequences of discrete solutions: for each n > 1, ex¬ 
tend u n by bilinear interpolation (or any other local interpolation procedure) on Q n := {x E 
d(x, dO) > n~ l Vd}i and dehne u*,u* : —> M by 

u*(x) := limsup(x), u*(x) := liminf u n (x). (23) 

n—>oo n—>CXD 


Lemma 2.14. and u* are respectively a sub-solution and a super-solution of \19\ right). 
Proof. Inspection of the proof of Corollary |2.9| yields the quantitive estimate 


mm a — 

an 


— y 

nd 


p{x) 


n" oj d 

. X&X n 


diam(fl) < u n < maxu. 

on 


(24) 


Note that n d Ylxex n p(x) f n p as n —> oo. (Also, the multiplicative term h* appearing i 


m 


Corollary 2.9 equals 1 with our choice of discrete domain X n = Q n n 1 Z“, since \x,x + e/n] n 


(X n U dfl) is non-empty for all x G X n . e G Z rf .) 

The maps u n G U n are therefore bounded independently of n, hence u* and u* are well 
defined. From this point, the announced result follows by a standard argument: let x G 0 and 
tp G C 2 {Pl) be such that u* — tp attains a strict global maximum at x G fh For each n > 1, let x n 
be the element of X n which maximizes u n — p. By monotonicity T> n p(x n ) > VnUn^Xn) = p(x n ). 
Observing that x n — > x as n —> oo, we obtain D(X 2 ip(x)) > p(x) in the limit. This establishes 
that u* is a super-solution, and likewise it* is a sub-solution. □ 


Corollary 2.15. One has u* = u*, and this map is the unique solution of (19 right). As a 
result, u n converges pointwise to Uoo as n —>• oo. 

Proof. By construction u* < u*, but by the comparison principle u* > u*, see Proposition 


2.13 Therefore ( |19| right) admits the solution Uoo := u * = u* which, again by the comparison 

oo: lirri sup u n (x) = 

□ 


principle, is its unique solution. Finally, for any x G Pi on has as n 
u*(x) = u*(x) = liminf u n (x), hence Uoo(x) = limu n (x). 
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3 Numerical experiments 


We implemented the thrccj^] Monge Ampere discretizations described in the introduction: the 
Finite Differences (FD) scheme T> FD , the Wide Stencil (WS) scheme T>g S , and the proposed 
scheme T>y. All presented experiments are three dimensional. The latter two schemes require 
choosing a collection B of triplets of orthogonal vectors, or a stencil V, see Table[l] We recall that 
scheme FD is consistent but not monotone (or degenerate elliptic). Scheme WS is monotone and 
thus benefits from the associated convergence guarantees, but suffers from significant consistency 
errors, see Figure [2] The proposed scheme is simultaneously monotone and consistent, provided 
the PDE solution hessian condition number is bounded, and the scheme stencil V is sufficiently 
large, see Proposition |1.7| and Corollary |1.8[ A consistency error arises when these conditions 
are not satisfied, see Figure [2j 

We limit our attention to synthetic test cases, posed on the unit cube fl :=]0,1[ 3 . A known 
convex function U : P —> M is numerically recovered from its hessian determinant p := det(V 2 t/), 
and its boundary values a := U\qq. Three test cases are considered. 

• (Quadratic) U(x) := ^(x,Mx), where M was chosen randomly, with \/||M||||M -1 || ~ 8.5. 


(Smoothed cone) U(x) := \J b 2 + \\x — xq|| 2 , with S := 0.1 and xq := (1/2,1/2,1/2). 


• (Singular, |FQ13| 1 U{x) := — -\/3 — ||x|| 2 . 

A Damped Newton solver is applied to the discrete system ([2]), starting from the trivial seed 
u(x) := ||x|| 2 . 

Quadratic test case. The chosen matrix M does not belong to the consistency set of schemes 
T>v and lbg S , with the chosen V and B\ hence the numerical error reflects their consistency error, 
and is resolution independent. On the topic of consistency, it was determined using Proposition 


1.7 and some semi-definite programming that the consistency set of the proposed scheme Vy 
contains all M E such that TY(M)/(detM)3 is less than 7.8 with the small stencil V, and 
less than 11.9 with the large V, see Table [T] The always consistent but non-monotone scheme 
V FD finds the exact solution for a range of resolutions, up to machine precision, but switches to 
a completely erroneous solution at other resolutions. 

Smoothed Cone test case. The test function U is smooth, but the test is not as simple as 
it looks since (i) VI/ varies quickly close to the center xq, but (ii) V 2 // is almost degenerate 
far from xo. Point (i) favors small stencils, while point (ii) requires a good angular resolution. 
Scheme WS is most efficient with a small triplets collection B at low resolutions, but with a 
medium or large one at high resolutions. Picking a stencil V is easier with the proposed scheme: 
the larger one here always yields a smaller error, albeit at a higher computational cost. The 
non-monotone scheme FD passes this test as well. We give some computation times for this test, 
which are typical of our experience. On a 2.7 Ghz Laptop using a single core, using a 50 x 50 x 50 
discretization grid X, computations took (in minutes): 4.7 and 9.7 with the proposed scheme, 
small and large V respectively; 6.4, 24 and 65 with scheme WS, small medium and large B 
respectively; 6.5 with scheme FD. 

Singular test case. The function U is not smooth at the point (1,1,1), where its gradient is 
formally (+oo,+oo,+oo). The non-monotone scheme FD completely fails this test: numerical 
error does not decrease as resolution increases, and initializing the Newton solver with the exact 
solution U does not even help, as observed in |FQ13j . Scheme WS will recover the solution, if 
the bases collection B is large enough, but numerical error decays quite slowly. The proposed 


lr The filtered method [FQ13 I. which accuracy may be competitive, was not implemented due to lack of time. 
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Scheme WS 

Orthogonal triplets #(^) 

Small 

Medium 

Large 

All within { — 1, 0,1 } 3 6 

All within {-2,- •• ,2} 3 43 

All within {—3, • • • , 3} 3 82 


Proposed 

Stencil vectors #(V) 

Small 

Large 

(1,0,0),(1,1,0),(1,1,1) 13 

Same and (2,1,0), (2,1,1) 37 


Table 1: Left: Stencil V used with the proposed scheme T>y. In addition to those indicated 
(oq, « 2 , 03 ), the stencil contains all vectors (eia^), £ 20 , 7 ( 2 ) i e 3 a o-( 3 )) obtained by permuting and 
changing the sign of their coordinates. Opposite vectors are identified when evaluating jf(V). 
Right: Collection B of orthogonal triplets used with the Wide-Stencil scheme Dg S . Triplets 
are counted up to reordering their components, and changing their sign: (eo,ei,e 2 ) ~ 
(eie CT (i), S 2 &a( 2 )i £ 3 e <r( 3 ))• Vectors with non-coprime coordinates, such as (2,2,0), are rejected. 

scheme works flawlessly: it has the same convergence guarantees as scheme WS, but yields an 
L°° error about 20 times smaller on a 50 3 grid. Curiously, numerical error is identical with the 
small and the large stencil V. 
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Figure 2: Relative consistency error ( T>um ~ det (M)) /T>um (blue = 0, red = 1) for a quadratic 
(|8]) map um associated to a matrix M = M{v), with v in the unit sphere (hence the spherical 
plot). First line: M(y) has eigenvalues 6 2 ,1,1, the former with eigenvector v. Second line 
likewise with eigenvalues 6~ 2 ,1,1. Third line M{y ) = R{y)DR(v) where R(v) is the rotation of 
axis v and angle 7r, and D is diagonal of entries 6,1,1/6. 



Figure 3: Log-Log plot of the L°° error, as a function of resolution n, for the three test cases. 
Discretization set ICO has n 3 points. Scheme stencils on Figure [l] 
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